##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(ggplot2)
library(parallel)

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--igc_cluster_file"), type = "character"),
    make_option(c("--dgc_cluster_file"), type = "character"),
    make_option(c("--diff_gene_file"), type = "character"),
    make_option(c("--pathway_path"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    sample_list_file <- "~/20220915_gastric_multiple/rna_combine/analysis/config/tumor_normal.list"
    igc_cluster_file <- "~/20220915_gastric_multiple/rna_combine/analysis/images/mfuzz_v2/mfuzz_plot_IGC.annotation.tsv"
    dgc_cluster_file <- "~/20220915_gastric_multiple/rna_combine/analysis/images/mfuzz_v2/mfuzz_plot_DGC.annotation.tsv"
    diff_gene_file <- "~/20220915_gastric_multiple/rna_combine/analysis/images/DiffGene/DiffGene.HugoSymbol.tsv"
    pathway_path <- "~/ref/Pathway/"
    out_path <- "~/20220915_gastric_multiple/rna_combine/analysis/images/mfuzz_v2"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
diff_gene_file <- opt$diff_gene_file
igc_cluster_file <- opt$igc_cluster_file
dgc_cluster_file <- opt$dgc_cluster_file
out_path <- opt$out_path
pathway_path <- opt$pathway_path

dir.create(out_path , recursive = T)

##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_igc_cluster <- data.frame(fread(igc_cluster_file))
dat_dgc_cluster <- data.frame(fread(dgc_cluster_file))
dat_diff <- data.frame(fread(diff_gene_file))

##########################################################################################

fisher_result <- function(P0.05_all,bg_gene,gene_pathway,pathway_name){
    ta <- length(which(P0.05_all%in%gene_pathway))
    tc <- length(which(bg_gene%in%gene_pathway))
    tb <- length(P0.05_all) - ta
    td <- length(bg_gene) - tc
    data_fisher <- matrix(c(ta,tb,tc,td),nrow=2)
    fisher_res <- fisher.test(data_fisher)
    
    gene_P_in_pathway <- P0.05_all[which(P0.05_all%in%gene_pathway)]

    if(length(gene_P_in_pathway)==0){
        gene_paste <- ""
    }else if(length(gene_P_in_pathway)==1){
        gene_paste <- gene_P_in_pathway
    }else{
        gene_paste <- gene_P_in_pathway[1]
        for(i in 2:length(gene_P_in_pathway)){
            gene_paste <- paste0(gene_paste,",",gene_P_in_pathway[i])
        }
    }

    res <- data.frame(pathway=pathway_name,enrich_top=paste(ta,length(P0.05_all),sep="|"),enrich_all=paste(tc,length(bg_gene),sep="|"),bggene_in_pathway=tc,pathway_gene=length(gene_pathway),P=fisher_res$p,OR=fisher_res$estimate,gene=gene_paste)
    res
}

computPathway <- function(pathway_path = pathway_path , P0.05_all = P0.05_all , bg_gene = bg_gene , cls = cls , output = output ){
    
    dir.create(output , recursive = T)

    for( pathway_file in grep( "symbols.txt" , list.files(pathway_path) , value = T) ){

        file <- paste0(pathway_path , "/" , pathway_file)
        dat_pathway <- fread(file)

        pathway <- unique(dat_pathway$pathway)
        pathway_res <- data.frame(rbindlist(mclapply(pathway,function(x){
            # print(x)
            gene_pathway <- subset(dat_pathway,pathway==x)$gene
            #gene_pathway_procod <- gene_pathway[which(gene_pathway%in%bg_gene)]
        
            res_enrichment <- fisher_result(P0.05_all,bg_gene,gene_pathway,x)
            return(res_enrichment)
        },mc.cores=20)))
        results <- subset(pathway_res,OR>1)
        resultso <- results[order(results$P),]
        resultso$P_fdr <- p.adjust(resultso$P,method="fdr",n=nrow(resultso))

        out_file <- gsub("[.]symbols[.]txt" , "" , pathway_file)
        write.csv(resultso,paste0(output,"/",out_file,".cluster",cls,".csv"),row.names=FALSE)
    }
}

diffCluster <- function(cluster_list = cluster_list , dat_use = dat_use , tumor = tumor){
    for( cls in cluster_list ){
        print(cls)

        use_dat <- subset( dat_use , CLUSTER == cls )

        P0.05_all <- unique(use_dat$Hugo_Symbol)
        bg_gene <- unique(dat_diff$Hugo_Symbol)

        output <- paste0(out_path , "/" , tumor , "_cluster"  )

        computPathway(pathway_path = pathway_path , P0.05_all = P0.05_all , bg_gene = bg_gene , cls = cls , output = output )
    }
}

##########################################################################################
## 不同的cluster分别聚类

dat_use <- dat_igc_cluster
tumor <- "IGC"
cluster_list <- unique(dat_use$CLUSTER)
diffCluster(cluster_list = cluster_list , dat_use = dat_use , tumor = tumor)

dat_use <- dat_dgc_cluster
tumor <- "DGC"
cluster_list <- unique(dat_use$CLUSTER)
diffCluster(cluster_list = cluster_list , dat_use = dat_use , tumor = tumor)


##########################################################################################
## 关注HallMarks通路
dat_plot <- c()
for( tumor in c("IGC" , "DGC") ){
    print(tumor)

    input_path <- paste0(out_path , "/" , tumor , "_cluster"  )

    for( cls in cluster_list ){
        file <- paste0( input_path , "/h.all.v7.5.1.cluster",cls,".csv")
        dat_tmp <- data.frame(fread(file))
        dat_tmp$CLUSTER <- paste0( "cluster" , cls)
        dat_tmp$Tumor <- tumor
        dat_plot <- rbind(dat_plot , dat_tmp)
    }
}

## 显著通路
dat_plot <- subset( dat_plot , P_fdr < 0.05 )
dat_plot$pathway <- gsub( "HALLMARK_" , "" , dat_plot$pathway )
dat_plot$CLUSTER <- factor( dat_plot$CLUSTER , levels = paste0( "cluster" , cluster_list[order(cluster_list)]) )

pathway_order <- unique(dat_plot[order(dat_plot$P_fdr , decreasing=T),"pathway"])
image_name <- paste0( out_path , "/Mfuzz.Hallmarks.pdf" )
dat_plot$pathway <- factor( dat_plot$pathway , levels = pathway_order , order = T)

p <- ggplot(dat_plot, aes(x=CLUSTER, y=pathway)) + 
  geom_point(aes(size=-1*log10(P_fdr),color=OR))+
  scale_colour_gradient(low="blue", high="red")+ 
  theme_bw() +  #设置背景
  facet_grid(.~Tumor) +
  labs(
    color=expression(OR),
    size="-log10(Q.value)") +
  theme(panel.background = element_blank(),#设置背影为白色#清除网格线
    legend.position ='right',
    panel.grid.major=element_line(colour=NA),
    plot.title = element_text(size = 8,color="black",face='bold'),
    legend.text = element_text(size = 10,color="black",face='bold'),
    axis.text.y = element_text(size = 10,color="black",face='bold'),
    axis.title.x = element_text(size = 10,color="black",face='bold'),
    axis.title.y = element_text(size = 10,color="black",face='bold'),
    axis.ticks.x = element_blank(),
    axis.text.x = element_text(size = 10,color="black",face='bold') ,
    axis.line = element_line(size = 0.5)) 
ggsave( image_name , p , height=10, width=12) 